Multi-patch methods in general relativistic astrophysics - I. Hydrodynamical flows on fixed 
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Many systems of interest in general relativistic astropiiysics, including neutron stars, accreting com- 
pact objects in X-ray binaries and active galactic nuclei, core collapse, and coUapsars, are assumed to 
be approximately spherically symmetric or axisymmetric. In Newtonian or fixed-background rela- 
tivistic approximations it is common practice to use spherical polar coordinates for computational 
grids; however, these coordinates have singularities and are difficult to use in fully relativistic models. 
We present, in this series of papers, a numerical technique which is able to use effectively spherical 
grids by employing multiple patches. We provide detailed instructions on how to implement such 
a scheme, and present a number of code tests for the fixed background case, including an accretion 
torus around a black hole. 

PACS numbers: 95.30.Lz, 47.ll.Df, 04.25.Dm 
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I. INTRODUCTION 

Numerical simulations have become one of our 
most valuable tools in building and refining mod- 
els of compact astrophysical objects and their envi- 
ronments, which are commonly associated with high- 
energy events like gamma-ray bursts, active galactic nu- 
clei and jets. X-ray binaries and supernovae. These 
rather exotic systems form one of the arguably most ex- 
citing physical laboratories we know of, where general 
relativity, nuclear physics, transport of radiation, mag- 
netohydrodynamics (and so on) contribute to their dy- 
namical properties. 

The rapid advance of computing performance has 
made it possible to simulate increasingly sophisti- 
cated problems. ^ But even with current high- 
performance supercomputers building general relativis- 
tic three-dimensional models poses a formidable chal- 
lenge, since the complexity of solving an hyperbolic 
problem of dimension n scales with the linear spatial 
resolution h like 0(l//i)"+^. It is therefore imperative 
to investigate the application of advanced, efficient nu- 
merical techniques in astrophysical models. 

This series of papers will focus on a particular ap- 
proach in which the computational domain is under- 
stood as a manifold covered by several distinct coordi- 
nate maps called patches. Above all else, this approach 
admits to cover spheres with smooth, and in particular 
singularity-free, coordinates, and also allows to employ 
many of the advantages of spherical polar grids, like 



^ For some recent results in general relativistic flow models see (ll,|^ 
lllllIIE[l,i,|l,[I^[Ill[ll[ll, though this list is far from being 
complete. 



the decoupling of radial and angular resolution and in- 
trinsically spherical domain boundaries, without shar- 
ing their major disadvantages. 

Multi-patch techniques have been used in different 



therein). The multi-patch irtfrastructure used for this 
paper ll27ll admits both overlapping patches and abut- 
ting ones (also called blocks), with the latter typically 
used when solving systems of equations with smooth 
solutions. In those cases the inter-block boundary in- 
formation is transported by the simultaneous approxima- 
tion technique | |28|] , combined with high-order finite dif- 
ference operators satisfying an algebraic property called 
summation by parts and associated dissipation operators. 
With these techniques the evolution system at the block- 
interfaces is decomposed into its characteristic struc- 
ture and any numerical mismatch between modes at 
these interfaces is dissipated away in time and/ or with 
increasing resolution in a numerically stable way (see 
m for more details). It has been possible to per- 
form very accurate simulations of scalar fields on a Kerr 
background jsill and fully non-linear simulations of dis- 
torted black holes | |32l| with this approach. 

As mentioned, these multi-patch techniques are de- 
signed for systems with smooth solutions, as is usually 
the case when solving the Einstein equations in vacuum. 
To model a larger class of astrophysical objects of inter- 
est one needs to deal with the presence of matter and 
the development of shocks in multi-patch simulations, 
which is precisely the goal of this series of papers. 

Here we treat hydrodynamical flows in the so-called 
test-fluid approximation, i.e., on a specified spacetime 
geometry. While this has interesting applications of its 
own, like the dynamics of non-self gravitating accretion 
disk models around black holes fill. l33l. Isl, ISa l and the 



normal mode spectrum of isolated neutron stars 13611 , we 
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will focus here on test cases with well-defined error func- 
tions (either Riemann problems or stationary solutions), 
to demonstrate the ability of the code to transport shock 
fronts across interfaces, and to keep stationary solutions 
near equilibrium. In later papers in the series we will 
consider the technique for the fully coupled evolution 
system and magnetized flows. 

This paper is organized as follows: In Section HH we 
give an overview of the multi-patch setup, the evolu- 
tion system, and the discrete techniques we are using. In 
Section Uni we give a specific description how to imple- 
ment such a scheme, which should be helpful for prac- 
titioners who would like to adopt this approach. In Sec- 
tion|IV]we discuss a number of multi-patch systems use- 
ful for test and production simulations. In Section IVl we 
present the results from a number of code tests, includ- 
ing shock tubes, rotating stars and an accretion torus 
around a black hole. Finally, in Section rvTl we summa- 
rize our results. 



II. THEORY AND DISCRETE TECHNIQUES 
A. The multi-patch setup in general relativity 




Figure 1: The general-relativistic multi-patch model used in 
this paper. The space-time sub-domain Mq C M of interest is 
assumed to be covered by a single global coordinate system cpQ. 
The target domain (Pg{Mq) C IR'* is then further covered by 
several patches associated with the charts (pci' which give rise 
to the local coordinate systems (p^. 



In general relativity (GR) one starts with a spacetime 
manifold (M, g), where M denotes the set of events and 
g the metric tensor field. Since the manifold is endowed 
with a differentiable structure D, an obvious choice for a 
continuum multi-patch model is, given an atlas [A] E D, 
a subset C A of charts : M — > ]R^ which cover the 
domain of interest (see, e.g., IIstII '). 

In our computational setup we will have, for practi- 
cal reasons, an additional element. We will assume that 
we are interested in solving our system of equations in 
a spacetime region Mq C M which can be covered by 
a single chart (pQ. That is, we will assume that there is 
a global coordinate system. This requirement is not fun- 
damental to using multi-patch techniques, but makes it 
easier to set up initial data, visualize results, and trans- 
port information between local patches, as discussed be- 
low. If needed, the assumption of such global system 
can be eliminated without any of the techniques of this 
paper changing. On the other hand, whenever we as- 
sume the existence of such global system, for definite- 
ness we will typically choose it to be of "Cartesian" type. 

In differential geometry language, a differentiable 
structure Di can be associated with the global region 
of interest ^q{Mo). A multi-patch setup, given an at- 
las [Ai] E Di, is then a subset of Ai covering (^g(Mo). 
For any <pQi E Ai, the product chart cpi = (pcL ° de- 
fines a local coordinate system on that local patch. Fig. [T] 
illustrates this relationship. 

A point of some technical relevance is the use of lo- 
cal or global tensor bases. At each point in Mq, tenso- 
rial quantities can be written in terms of global or local 
coordinates (in differential geometry language, the tan- 
gent space bases associated with (pQ or (pi^, respectively. 



can be used). The latter appears to be the more natural 
choice and will be employed to represent the hydrody- 
namical variables in this paper to enforce conservation 
when this is available. A global system is still useful 
(though not strictly necessary) to serve as a "reference" 
frame. On the other hand, there is no obvious conserva- 
tion law for the metric sector, and therefore there is no 
advantage in using local coordinates. In fact, it is eas- 
ier to use global ones, and therefore when we solve for 
the metric we do so for its global components (see, e.g., 
iH-HMl) • In such case, a standard finite-difference dis- 
cretization delivers partial derivatives consistent with 
(pi^ which are then transformed by the Jacobians asso- 
ciated with (pQi. 



Once the multi-patch system has ben set up, the sys- 
tem of equations to be solved can be discretized on 
each patch and two types of boundary conditions ap- 
pear: those associated with the outer boundaries of the 
whole domain Mq, and those connecting neighboring 
patches. In our applications each point in Mq will rn 
general be associated with the interior of exactly one 
patch, and with the boundary region of several (two or 
more) patches (see Fig. [JJ. Multi-patch techniques can 
generally make use of the overlap region to define ghost 
zones, or they use the interfaces to establish a common 
surface for the communication of variables. 
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Figure 2: The notions of interior, interface and boundary re- 
gion associated with a patch. The global coordinate domain 
i^c(Mo) is decomposed into the patch interiors as indicated in 
the figure, with interfaces adjoining two interiors. The charts 
(pQj^ therefore define an interior and a boundary region in 
cp'liMo). 



B. General relativistic hydrodynamics 

1. Evolution system 

We assume an energy-momentum tensor of the usual 
form Ii38i1 



{v{zv)) = s{v{zu)). 



(6) 



Here, the tupel v = [p, u, u^, uV, u^, P) denotes a choice 
of primitive variables, which, in general, are an implicit 
function of the conserved variables w. The fluxes F' (v) and 
sources s{v) are expressed in terms of the primitive vari- 
ables, which makes it necessary to obtain the primitive 
from the conserved variables at each evolution step. 

For calculating the primitive variables, we make use 
of the "2D scheme" from jioll , since it easily lends itself 
to the inclusion of magnetic fields. The conversion is a 
non-linear root-finding problem, which is being solved 
using a Newton-Raphson scheme, and as such is one 
of the most delicate parts of the implementation. The 
Newton-Raphson scheme only converges given a suffi- 
ciently close initial guess: We use the value of the primi- 
tive variables at the last sub-step for this purpose. Also, 
in some cases the conserved variables could potentially 
obtain values which are not compatible with any set of 
physical primitive variables and, finally, the round-off 
errors may lead to subtle problems, e.g. a Lorentz factor 
of w 1 — 10~^^. Most of this is discussed in lioll and in 
the accompanying source code examples. 

The pressure function P = P{p,u) is chosen according 
to the physical properties of the fluid, which is in gen- 
eral determined by the equation of state. For purposes 
of this publication, we will assume the pressure to be 
obtained from the gamma law P = (F — 1)m. 



T"^ = {^p + u + P)u''u^ + Pg''^ (1) 

where p is the fluid's rest-frame mass density, u is 
the internal energy density, li" the four-velocity, P the 
isotropic pressure and g"^ the contravariant components 
of the spacetime metric. Given a patch with coordinate 
system (pi the conservation laws of mass and energy- 
momentum 

Vaipu') = (2) 
VaT^ = (3) 

are transformed into an evolution system by a 3+1 split 
(see, e.g., jH]) 

MV^Pu')+M^/^P^n = (4) 

HV^8^'a)+MV^gTa) = V^gT'd^'ac (5) 

where f is the time coordinate, i is a local space coor- 
dinate, g = det{g„i,) is the determinant of the space- 
time metric, and V^ac are the Christoffel symbols asso- 
ciated with it. We therefore promote the expressions 
D = \J—gpu* and Qa = \/—gT*a to evolution variables, 
and collect them into the tuple iv = (D, Qa). The equa- 
tions can then be expressed in the flux conservative form 



2. Discretization 

On any patch, we use a standard finite-volume 
scheme 14111 to update the evolution variables. The com- 
putational domain is broken into blocks, and each grid 
point is assumed to coincide with the center of a finite 
cell. The primitive variables associated with a grid point 
are naturally interpreted in terms of volume averages 
over the cell. The evolution system eqn. [6] can then be 
integrated over the cell to yield the weak form of the 
equations, which represents the update of the volume 
averages in terms of fluxes across the cell interfaces and 
source terms: 

dt / wdadbdc (7) 

i(fll,bl,ci) 

= - L ^ (FHv)\a,-FHv)U dbdc 

f {12,1:2) / r, T \ 

- / (f'{v)\,^-FHv)U dadc 

- [f\v)\c,-F\v)\cA dadb 

r(a2,h,C2) 
+ I s{v) dadbdc 
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The fluxes across the cell interfaces are obtained in 
two steps: First, the primitive variables are extrapolated 
to their (left and right) interface values using a recon- 
struction algorithm, and then the thus-defined local Rie- 
mann problem is solved with an approximate method. 
For purposes of reconstruction, we use the MC (mono- 
tonized central) algorithm. If m, is the value of the prim- 
itive variables at zone i, we define the reconstructed 
quantities vi and vr at the interface location / + ^ as ll42ll 



VL = Vi + Ai/2 

vr = - A,-+i/2 

sgn(A,-+i)min(2|A,-|,2|A,-+i|, 
|A,- + A,-+i|/2) AA-+i>0 
otherwise 



(8) 



A/ = 



The resulting local Riemann problem defined by the 
states Vi and Vr is approximated using the HLL (Harten, 
Lax, van Leer) flux formula jisll 



P = {c,ni„P{vR) + C^a^F'{vL) (9) 
{wr -Wi)){ 

^min^niax) 

Cmin = -min(0, c^,c^) 
Cmax = max(0,cj,cj) 




Figure 3: Illustration of the patch interface boundary treat- 
ment. The interface, here represented by the vertical thick blue 
line in the center, separates patch L {left) and patch R {right). 
The coordinate lines of both patches are represented in the lo- 
cal coordinates of patch R, which makes the lines of L appear 
curved in general. To evolve the system on patch R, bound- 
ary gliost zones, here represented by non-filled circles, are intro- 
duced as a linear extrapolation of patch R's coordinate lines, 
and, before the time update is performed, receive data from a 
data interpolation operation on patch L. Afterwards, the data 
is transformed to the tensor basis defined by patch R's coordi- 
nate system via the transformation map. 



Here and c+ are the minimal and maximal char- 
acteristic speeds associated with the variables Vi r, and 
zvi^R are the conserved variables obtained from Vi^r. 
This flux formula has the advantage that it does not re- 
quire the full characteristic decomposition of eqn. |5j and 
is thus well-suited to be extended to the more compli- 
cated case involving magnetic fields Q. The charac- 
teristic speeds are obtained from the rest-frame sound 
speed of the fluid by a Lorentz transformation to the lo- 
cal coordinate frame jisll . 

A particular difficulty in solving general relativistic 
hydr ody namics problems are regions of very low den- 
sity ll46l I47I1 , e.g. those outside a neutron star. In those 
regions, we use an artificial atmosphere of low density 
to make the problem tractable. This atmosphere effec- 
tively acts as a boundary condition on the stellar mate- 
rial. 



3. Boundary treatment 

To use an xmmodified finite-volume scheme also at 
the boundaries of the interior patch domain, we fol- 
low the same approach typically used in adaptive mesh 
refinement implementations: we introduce ghost zones 
in the boundary region which are updated using an 
interpolation scheme. The results in this paper have 
been obtained with a first-order operator, which does 



not introduce new extrema and is compatible to the 
(at most) second-order accurate MC reconstruction tech- 
nique. The patch interface boundary treatment is illus- 
trated in Fig. [3l 

Since the patches have different local coordinate sys- 
tems, we need to use the transformation maps between 
the patches to obtain the local representation of inter- 
polated quantities. Internally, first a transformation to 
the global coordinate system is performed, followed by 
a separate transformation to the appropriate local sys- 
tem. For the case of general relativistic hydrodynam- 
ics on specified backgroimds, we interpolate the set of 
primitive variables for purposes of interface reconstruc- 
tion. These quantities are components of tensors and 
therefore subject to trivial transformation laws IstIi . 

The decomposition of each patch into sub-domains, 
for purposes of implementing a distributed comput- 
ing model, introduces additional boundaries which are 
treated by introducing ghost zones. No interpolation 
or coordinate transformation is needed in this case, and 
the synchronization operation copies boimdary data be- 
tween these sub-domains. 

The outer boundaries of the computational domain 
are also handled with ghost zones. For purposes of im- 
posing an outflow boundary condition, we copy data 
from the first cell inside the computational domain to 
the ghost zones [48], though we will typically prefer to 
select a grid setup where material ejected from it is en- 
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tirely contained on the grid during the course of the evo- 
lution. In this case, the artificial atmosphere is the effec- 
tive boundary for the fluid. If material should leave the 
computational domain, the exact form of the boundary 
condition is important only in those cases where there is 
a significant back-reaction of the material on the central 
object, or other regions of interest. 

4. Setup of initial data 

For convenience, initial data is always expressed first 
in terms of the global coordinate system, and mapped 
to the patch locations in this basis. Afterwards, the co- 
ordinate transformation to the appropriate local basis is 
performed if the quantity is not a scalar field. All ini- 
tial data specifications in this paper will therefore also 
be written in an appropriate global system. 

C. Computational framework 

The code has been implemented as a module into the 
Cactus computational framework |49l,l53| . and using the 
Carpet driver jsil l52ll to administrate the multi-patch 
infrastructure. In our implementation (as stated above), 
the Carpet driver reserves a number of boundary ghost 
zones beyond each patch boundary, which are then filled 
using interpolation and tensor transformation. For a 
first order interpolation scheme, each boundary ghost 
zone is located in exactly one cell in the interior of an- 
other patch. However, if the scheme was to be extended 
to higher order interpolation operators, a non-central 
stencil would need to be applied. 

D. Limitations of our approach 

In this section, we collect the most important known 
limitations of the approach followed in this paper. 

The most obvious limitations are clearly those re- 
lated to the physical model: the assumption of a fixed 
spacetime neglects motions of the gravitational field and 
therefore suppresses all gravitational instabilities, and is 
only strictly justified in code tests involving equilibrium 
systems. Our fluid model, besides making the assump- 
tions common for all hydrodynamical approximations, 
models one fluid component without tangential stresses, 
and our choice of the pressure relation neglects micro- 
physical properties. This paper does not include mag- 
netic fields either, although we will turn to full GRMHD 
in a subsequent publication. 

On the numerical side, we have a number of restric- 
tions. While many of these are common to many ap- 
proaches in the field, we find it useful to give a detailed 
list here: 

• The scheme that we use is at best second order ac- 
curate. In particular, it drops to local first order 



near maxima and stellar surfaces. The global con- 
vergence order might be less than two in practi- 
cal applications. This approach is very robust, but 
may be too inaccurate for very long term simula- 
tions (hundreds of dynamical timescales), or tur- 
bulent phenomena like the MRI 115311 . There are 
higher order schemes available IdHI , and we are 
evaluating using such techniques.^ 

• The inter-patch boundary treatment does not triv- 
ially lend itself to higher order implementations, 
since we need a non-oscillatory higher order inter- 
polation operator to fill the boundary ghost zones. 
A compromise would be to use Lagrangian opera- 
tors, but drop to trilinear interpolation if unphysi- 
cal variables are produced, or use an ENO scheme. 

• The outer boundary conditions are only approxi- 
mate, and may lead to unphysical variables (see 
l39ll for an improved scheme with densitized vari- 
ables). This may have practical and fundamental 
consequences. On a practical level, undesirable ar- 
tifacts may appear at the boundaries, though in 
stationary systems, those are only related to arte- 
factual flows leaving the compact support of the 
equilibrium object. However, if the system does 
show physical outflows across the boundary, those 
will be modeled only approximately. For super- 
sonic or even causally disconnected flows, how- 
ever, this restriction may not be dynamically rele- 
vant 

• The artificial atmosphere is used to employ a uni- 
fied computational scheme. However, to produce 
an exact discrete representation of an equilibrium 
star or disk, the discrete scheme would need to be 
modified near the surface. The addition of atmo- 
spheric fluid usually introduces error levels which 
are significantly smaller than errors from the inter- 
nal dynamics of the star in real applications, but 
for equilibrium problems, which may be domi- 
nated by the surface errors, the code may not show 
convergence to the correct solution even when re- 
ducing the atmospheric density. 

• The transformation from conserved to primitive 
variables works well in most cases, but can be a 
cause of considerable difficulty in some situations. 
This problem may be more severe when using a 
tabulated equation of state and coupled magnetic 
fields; in fact, the technique we currently use may 



^ Higher order reconstruction operators like PPM ts^ do not increase 
the overall accuracy but may decrease local error levels. We do have 
a simplified implementation of PPM in our code (without flatten- 
ing), but we did not find it superior in the particular test cases dis- 
cussed here. 
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need a considerable extension to also operate in 
these cases ll56l] . 

• The multi-patch setup, for all its advantages men- 
tioned above, requires a somewhat sophisticated 
software infrastructure, and a certain reduction 
in computational efficiency is introduced by the 
boimdary commimication. However, a multi- 
patch setup scales more favourably when com- 
pared to Cartesian mesh refinement when one 
considers increases in radial resolution or an in- 
crease of the computational domain. The com- 
putational cost of a multi-patch simulation scales 
as 0{N^) and 0{N), respectively, in these two 
cases, while a mesh refinement simulation scales 
as O(N^) and 0{N^), respectively The difference 
comes from the fact that the radial resolution re- 
mains unchanged in these cases in a multi-patch 
setup, which is not possible in a Cartesian mesh- 
refinement simulation. 

• The patch interface grid points from both sides 
must match, i.e., the grid topology introduces con- 
straints on the possible choices of patch cell num- 
bers. While this could be relaxed for the hydro- 
dynamical scheme alone, the SBP/penalty tech- 
niques used in the generalized harmonic code for 
solving Einstein's equations | [2^ require these con- 
straints, and we are ultimately interested in the 
full evolution system. 

III. IMPLEMENTING A MULTI-PATCH SCHEME 

In this section we describe how to implement a multi- 
patch scheme, either from scratch or by extending an ex- 
isting imigrid code. Multi-patch techniques can be use- 
ful for any numerical code which evolves systems with a 
certain symmetry (for example, single stars or stars with 
accretion disks), mostly because the grid is adapted to 
the problem, and the angular and radial resolutions are 
decoupled. Because of these and other advantages, this 
section should be of interest to practitioners in the field 
of computational astrophysics, whether Newtonian or 
relativistic. 

What amount of work can one expect, and what ben- 
efits can be gained in practical terms? We will answer 
these questions in comparison to the two main alterna- 
tives to multi-patch grids: 

• In comparison to spherical polar grids, the multi- 
patch technique avoids the need to exclude the 
axis of symmetry and impose artificial boundary 
conditions there. This is important, in particu- 
lar, for systems where outflows are generated and 
collimated on or near the axis. Also, the typical 
finite-difference or finite-volume Courant limita- 
tions to the time step near the poles of each sur- 
face r — const are avoided. Three capabilities need 



to be added: geometric terms which are not hard- 
coded into the equations, several grids instead of 
one, and a linear interpolation and transformation 
at the boundaries (see below). The finite volume 
scheme, Riemann solver, and local physics are im- 
affected. 

• In comparison to mesh refinement, the multi-patch 
technique offers decoupled radial and angular res- 
olution, which is of particular use far away from 
the central object (e.g. to extract gravitational ra- 
diation at large distances). However, mesh refine- 
ment is superior for processes which do not have 
explicit approximate symmetries, e.g. violent in- 
stabilities or binary mergers. A mesh refinement 
code can already handle several grids and bound- 
ary interpolation, so the only capabilities needed 
in addition are the coordinate transformations at 
the boundary and the geometric terms for the local 
coordinates. The finite volume scheme, Riemann 
solver, and local physics are again unaffected. 

Given a certain finite volume / finite difference code, 
these steps need to be performed: 

• The code needs to be able to handle several sepa- 
rate grids. These grids are all logically Cartesian 
and independent, so in many cases this might just 
be an additional loop statement. The local evo- 
lution scheme (e.g. flux calculation, finite differ- 
ences, update) is logically unaffected. 

• Each grid cell needs to store (or calculate during 
run-time) its location in global coordinates, typi- 
cally Cartesian x' = (x, y, z), and in local coordi- 
nates, called a' = {a, b, c) below. In addition, the 
Jacobian }' j = / daJ and its inverse are needed; 
these can be obtained from the equations in Sec- 
tionHVl 

• The patches need to use local coordinates in a 
certain range; say [—1,+!] for the interior, with 
a' = ±1 on the boundary, and a number of grid 
points extruding beyond the boundary for setting 
interpolation data. This depends on the stencil; for 
the monotonized central scheme used here, we use 
two ghost cells. 

• Before each time update, and after the conversion 
from conserved to primitive variables, the bound- 
ary ghost cells need to be set by a (for example, tri- 
linear) interpolation operation. For the patch sys- 
tems presented in Section HVl it is known, for each 
patch, which other patch needs to be interpolated 
on,'^ so a simple list of six entries for each patch 



^ In general, this is not the case. Given the location of the ghost point 
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is enough to store the information. Given the tar- 
get patch, we transform the location of the ghost 
point into local coordinates on the target patch using 
the equations in SectionlTVl and use traditional tri- 
linear interpolation to get the data. 

• Although boundary data has now been set, it is 
still expressed in the coordinate system of the tar- 
get patch. Say, for relativistic hydrodynamics, we 
are interpolating boundary data for patch 4, and 
the data is coming from patch 5. After the in- 
terpolation, patch 4 has one side of its boundary 
stencil (the side adjoining to patch 5) filled with 
interpolated values for (p, m, m^^^, m^^^, u^^p. As 

indicated, the 3-velocity components are still ex- 
pressed in the local coordinate system of patch 
5. Therefore, a simple transition map m|^j — 

(3fl|^j/8x-')(3x-'79fl^5pM(5^ is applied at each ghost 
cell, where (dxi /da^^^^) is the Jacobian local 
global to transform m'^^j to global coordinates, and 
(8a|^j /8x-') is the Jacobian global — > local on patch 
4. 



M processes, in order to reduce commtmication 
overhead.^ 

IV. MULTI-PATCH SYSTEMS 

In this section we specify the particular set of multi- 
patch setups used in this paper. A general introduction 
to the kind of systems useful for single stars or black 
holes can also be found in ll29l,l30ll . 

For convenience and clarity, we will denote those 3- 
coordinates associated with the global coordinate sys- 
tem as x',i = 1...3 or (x, y, z), and those associated 
with a particular local coordinate system as a',i = 
1 ... 3. If necessary, we indicate the particular patch p 
with a subscript, as in a'^. 

A. The uni-patch system 

We refer to a setup with one patch and a coordinate 
transformation with a constant diagonal Jacobian of the 
form 



• The numerical scheme needs to be able to work on 
a general background metric. For relativistic codes 
that is already contained in the covariant form, but 
a Newtonian code requires an addition of the geo- 
metrical terms. 

• It is convenient to describe the initial data in terms 
of the global Cartesian coordinate system. That 
is, at each grid point all tensorial data can be 
represented in global coordinates {x,y,z) and af- 
terwards transformed to the patch-local coordi- 
nate systems using the standard tensor transfor- 
mation laws. In general relativistic hydrodynam- 
ics on fixed backgrounds, we need to transform 

(G) 

the 3-velocity W^q-^ and the 4-metric g^u from 
global into local coordinates on patch m: = 

{da^^^/dxi^^y and gjf) 



(m)'i>kl 



To parallelize the code, a simple domain decom- 
position technique can be used for each patch. In 
addition, it may be useful to have each patch dis- 
tributed to as few processes as possible. For exam- 
ple, for a setup with N patches one might want to 
use N X M processes and assign each patch onto 



ji,^dx^^^gi. 

'^-dai ' 



(10) 



as a uni-patch system, where j is the Kronecker sym- 
bol and c G ]R a number. The associated transformation 
between local and global coordinates is the affine map 



ca 



V,i 



(11) 



and therefore we only need to specify the scale c and the 
3-vector location of the origin V . 

This patch system most closely matches the unigrid 
setups used in some three-dimensional simulations, and 
since it contains no patch interfaces it is used to test the 
performance of the code in the patch interior.^ 



B. The distorted two patches system 

The distorted two patches system consists of two 
patches with coordinates 



Patch 0: x' = a\, f = 1 . . .3 



(12) 



in global coordinates, a generic approach is then to transform the 
point to all local coordinate systems, and check whether it is in the 
appropriate interior range, say [—1, +1]^. We perform this step once 
during initialization. 



^ Our load distribution scheme is in fact more complex and can also 
efficiently handle patches of different sizes. 

^ It is also easily possible to make this system toroidal by performing 
a topological identification between some boundaries, though we 
will not make use of this option in the examples below. 
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I 1 -0.2 
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Figure 4: An illustration of the distorted two patches system. 
The left patch is constructed using an identity map between 
local and global coordinates, whereas the right patch follows 
from eqn. [TS] The common interface is located ai x = 1, 
which also coincides with the local coordinates a\ = \ and 
a\ = \. Note the ghost zones extending beyond the interface 
(cf. Fig.H). 

Patch 1: xi = -^(4 + if + 3, ; £ R (13) 

x' = flj/ = 1 ■ • . 3, / 7^ /. 

The internal coordinate range of each patch is chosen 
to be (0, 1)'', and the common interface is then located at 

= fl^ = 4 = 1- The second patch is quadratically dis- 
torted and matches the first one at the upper boundary 
of coordinate j. An example with = 1 is illustrated in 
Fig. II 



introduces coordinate singularities. Therefore, either 
one opts to cover the sphere using at least to patches 
with properly rotated spherical polar coordinate sys- 
tems {yin-yang grid (S^), or one uses a patch system 
which does not derive from spherical polar coordinates. 
We will use the latter approach here. 

The cubed sphere patch system on a sphere is concep- 
tually constructed by imagining a cube con-central with 
the sphere, introducing the canonical Cartesian restric- 
tions on each surface of the cube, and projecting these 
coordinate lines onto the sphere. When mapping the ra- 
dial coordinate via the identity, this produces six patches 
with coordinate transformations 12911 



Patch 0: 


X = 


Go 

CO 


Go^o ^ _ GoflJ 
Eo ' Eo 


Patch 1: 


X = 


El ' 


Gi 

y = z = 


Gia\ 
El 


Patch 2: 


X = 


-G2 

E2'^ 




G2a\ 
E2 


Patch 3: 


X = 


T'' y 




G34 
E3 


Patch 4: 


X = 


-Gia\ 
£4 ' 




G4 
E4 


Patch 5: 


X = 


E5 


Gsflr 
£5 


-G5 
E5 








' + 






G,= 


4 [^0(1 










[-1.1] 







C. The cubed-sphere six patches system 

Whereas the two multi-patch systems discussed 
above are useful mostly for purposes of code testing, the 
cubed-sphere six patches system [291] is a setup which can 
be directly used to efficiently model black holes and ac- 
cretion disks aroimd black holes. The system covers a 
spherical region with a central sphere cut out {excised). 
The central excised sphere can be used to exclude part 
of the interior of a black hole, since it is causally discon- 
nected from the exterior region. In a fully relativistic set- 
ting, a symmetric h5^erbolic formulation of Einstein's 
equations with only physical speeds of propagation ad- 
mits to use the interior spherical surface as an outflow 
boundary, since the continuum property of causal dis- 
connection is then also represented on the discrete level. 
We will make use of this fact in the coupled evolutions 
in a later paper in this series. 

The domain is covered by a family of spheres with dif- 
ferent "radial" coordinates r = ^ simple 
approach would be to cover almost the entire sphere by 
spherical polar coordinates, which, imfortunately, also 



Here, and r\ are free parameters which determine 
the location of the two outer boundary spheres in terms 
of global coordinates. Patches to 3 are constructed in 
the equatorial plane counterclockwise starting from the 
positive X axis, and patches 4 and 5 intersect with the z 
axis. An illustration of this patch system in the equato- 
rial plane z = is given in Fig. [S] 



D. The cubed-sphere seven patches system 

While the six patches system is well-suited to excise 
the interior of black holes from the computational do- 
main, there are also situations (for example, stellar sim- 
ulations) where we would like to cover the entire inte- 
rior of a sphere with regular coordinates. One approach 
to achieve this is to cover the central part of the compu- 
tational domain with Cartesian coordinates (i.e., a uni- 
patch as described above), describe the outer boundary 
as a sphere, and decompose the region between the sur- 
face of the Cartesian patch and the outer boundary into 
six additional patches Isoll . Then, on the outer boundary. 



-2.5 



-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 



-2.5 



-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 



Figure 5: An illustration of the cubed-sphere six patches sys- 
tem. The diagram shows a cut of grid lines with the equatorial 
plane z = 0. In the notation of eqn. [151 patch is to the right 
(in red), counterclockwise continuing with patch 1 (in green), 
patch 2 (in blue) and patch 3 (in violet). The inner spherical 
boundary can be used as an excision surface for modeling black 
holes. 



spherical coordinates are fixed as in the six patches sys- 
tem, whereas the intermediate surfaces of constant coor- 
dinate (which are spheres in the six patches system) 
are deformed in a way to interpolate between a spheri- 
cal section and a flat face of the central cube. 

hi detail, the transformations between local and 
global coordinates which specify the cubed-sphere 
seven patches system are as follows (please note that we 
introduce the central cube as patch 6 to retain the count- 
ing convention used in the six patches system): 



Patch 0: x 

Patch 1: X 

Patch 2: x 

Patch 3: x 

Patch 4: x 

Patch 5: x 

Patch 6: x 



Go 



^0 



— Gia\ 

-G2 
F3 ' 



' y 



y 



Gi 

-02^2 

~F^' 
-G3 . 



mi y 



Fa ' 
_ Gjul 

Gsal 



ml, z 



Fo 

G\a\ 
G2a\ 

F3 
Gi 

' = Ti 
-G5 



(15) 



Figure 6: An illustration of the cubed-sphere seven patches 
system. The diagram shows a cut of grid lines with the equa- 
torial plane z = 0. In the notation of eqn. [I6l patch is to 
the right (in red), counterclockwise continuing with patch 1 
(in green), patch 2 (in blue) and patch 3 (in violet). The central, 
Cartesian patch is patch 6 (light blue). The parameters for the 
patch system were chosen to be rg = 0.5 and rj = 2, so the 
outer boundary is a sphere of radius 2 in this diagram. The 
outermost two grid lines, however, are ghost zones required 
for the MC reconstruction algorithm (see Section lTlB3t . 



Fi 



(ri-G,) + (G,-ro)£,V/' 



n - ro 

1^2 I i„2\2 



ml 



G, = i[ro(l-flf)+ri(l + flf)] 
a{e[-l,l] 

Here, tq and ri are free parameters which determine, 
in terms of global coordinates, the location of the outer 
boundary (r^) and the extent of the central cube (tq). An 
illustration of the cubed-sphere seven patch system in 
the equatorial plane z = is given in Fig. |6l 



E. The cubed-sphere thirteen patches system 

A (minor) disadvantage of the cubed-sphere seven 
patches system is that, beyond the Cartesian patch 
around the origin in global coordinates, all surfaces of 
constant fl? for i = ... 5 up to the outer boundary at 
fl? = 1 are not spherical, but rather "intermediate" be- 
tween spheres and cube surfaces. To cover a star with 
spherical surfaces, a different approach is to combine the 
advantages of the six and seven patches systems by cov- 
ering a large region of the star with a six patches system. 
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-4-3-2-101234 

Figure 7: An illustration of the cubed-sphere thirteen patches 
system. This system consist of an outer six patches system, 
where all surfaces of constant coordinate are spheres, and 
an inner seven patches system matched to the irmer boundary 
of the outer system to cover the origin in a regular manner. 



(Or = 1, Mr = 1.5 -10 ^, Mr = 

In all cases we assume an equation of state of the form 
P = (r — 1)m with r = 4/3, and the metric has the 
Cartesian Minkowski form g^,y = rj^,y The Sod test, in 
particular, will also be used for shock tube experiments 
on several patches below. 

For the uni-patch system, we choose the free parame- 
ter c in eqn.[10]to be unity, and we restrict attention to a 
one-dimensional problem where the initial contact sur- 
face is located at x = 0.5. The grids, therefore, are only 
specified by the number of cells in the x direction, where 
we choose values from 50 to 800. 

To generate a measure of error, we calculate the exact 
solution to the the special relativistic problems using the 
riemajin code \^^. The error function for each primitive 
variable is simply taken to be a norm over the difference 
function between the exact solution and the discrete re- 
sult. The boundaries are set to the exact solution pro- 
duced with riemainn. 

The results of evolution of the Sod, "Simple" and blast 
wave problem are illustrated in Figs. [3 19] and [lOl The 
code has some inherent dissipation due to the use of the 
MC limiter, but it converges to the exact solution. 



B. Sod test on the distorted two patches system 



and making the origin regular by introducing a seven 
patches system in a way that the parameters rj of the 
seven patches system and rg of the six patches system 
match. This construction can also be used, for example, 
to extract gravitational waves without needing to inter- 
polate to spheres in order to decompose the solution into 
its different multipole components. A cut of the coordi- 
nate lines with the equatorial plane z = is illustrated 
in Fig. [7] 



We have performed the Sod test (eqn. [T6ll on the two 
distorted two patches system described in Section HV Bl 
The first and second patch both have n x 1 x 1 cells, 
where n = 50 . . . 800, and the patch interface is located 
at X = 1 in terms of global coordinates. As before, the 
boundary ghost cells are set to the exact solution. 

Fig. [TT] demonstrates that the multi-patch approxima- 
tion leads to a convergent transmission of the shock 
across the interface, although both patches use different 
local tensor bases, and a non-trivial Jacobian is associ- 
ated with the second patch. 



V. RESULTS 
A. Shock tubes on the uni-patch system 

To investigate the code's ability to evolve special rela- 
tivistic Riemann problems without the additional com- 
plications of non-trivial Jacobians and inter-patch inter- 
faces, we perform three kinds of standard tests specified 
by the fluid states 



Sod test: 


PL = 


= 1, "L = 


1.5, 









PR -- 


= 0.125, Mr = 


0.15, 


"r = 


"Simple" test: 


Pi = 


-- 10, Ml = 


= 20, 


m1 = 


= 




PR -- 


= 1, Mr = 


10" 




= 


Blast wave test: 


PL = 


= 1, Ml = 


1.5- 


10^ 


"L = 



C. Sod test on the cubed sphere six patches system 

To perform the Sod test (eqn. [T6l l on the cubed sphere 
six patches system described in Section IIVCI we pre- 
pared a setup with n x n x n cells per patch, where 
n = 20,40,80, and chose the free parameters rg and ri, 
which specify the inner and outer boundary radius in 
terms of the global coordinates, to be given by rg = 1 
and ri = 2. The outer boundaries of the domain are set 
to the exact solution produced by riemann. 

The evolution of the density for the case n = 80 is 
shown in Fig. [12l This plot shows the density function 
in the equatorial plane, and the motion of the waves re- 
sulting from the Riemann problem. Global convergence 
to the exact solution, which is again constructed using 
riemann and then transformed to the local coordinate 
systems as required, is demonstrated in Fig. [l3l 
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Figure 8: Sod test on the uni-patch system. For each of the primitive variables p, u, and m^, the left plots show the comparison 
of the numerical result with the exact solution for different times (spaced in units of Af = 0.1), and right plots show the error 
function (/j norm of the difference between the exact solution and the discrete result) over time. 



As explained in Section III B 3[ the boundaries are 
treated by setting boundary ghost zones using an inter- 
polation operation before each time update, followed by 
a suitable coordinate transformation. Fig. [14] shows a 
pseudo-Schlieren plot of the transmission of the shock 
front across an interface. The code produces a number 



of minor reflections as expected, but the overall shape of 
the shock front remains intact. 
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Figure 9: "Simple" shock test on the uni-patch system. Same as in Fig.[8] 
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D. Sod test on the cubed sphere seven patches system 



The setup of this test is very similar to Section lVCI but 
uses the cubed sphere seven patches system described 
in Section [IV Dl The free parameters rg and ri, which 
specify the extent of the central cube and the location 
of the outer spherical boundary, are set to rg — 0.5 and 
ri = 2. Again, we use nx nx n cells, with n = 20, 40, 80, 



on each patch, and set the outer boundary ghost zones 
to the exact solution. 



As Fig.[T5]shows, the evolution proceeds very similar 
to the six patches system and is convergent to the exact 
solution. 
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Figure 10: Blast wave test on the urd-patch system. Same as in Fig.|8] 



E. Sod test on the cubed sphere thirteen patches system 

In the same way as in Sections IV CI and IV Dl we per- 
form the Sod test on the cubed sphere thirteen patches 
system described in Section flV El Since this setup con- 
sists of seven inner patches matched to six outer ones, 
we have three free parameters available. For purposes 



of this particular test, we use tq = 0.5, ri = 2 and r2 = 3. 
The results are presented in Fig. [I6l and show conver- 
gence to the exact solution. 
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Figure 11: Sod test on the distorted two patches system. This system, defined in Section HVBI consists of a left patch {patch 0) 
which is undistorted, and a right patch {patch 1) which has a non-trivial Jacobian associated with the transformation from local to 
global coordinates. The left panel shows the primitive variables p, u and i/^ in dependence on the global coordinate x, for different 
coordinate times {At = 0.2). The interface between the patches is located at x = 1. In the lower left diagram, note how the velocity 
component i/^ appears discontinuous across the interface since it is represented in different tensor bases on each patch. The right 
panel shows global convergence with respect to the li norm. (The left panel does not contain the graphs for n = 200 to enhance 
the visual clarity of the plots.) 
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Figure 12: Sod test on the six patches system: evolution of the density for the case where each of the patches has resolution 
80 X 80 X 80. The white lines indicate the boundaries of the six patches which are used to cover the computational domain. 
The yellow surface is the cut of the plane 2 = with the grid boundaries, and the density function is shown at z = 0. The left 
panel shows a perspective view of the evolution for coordinates times 0, 1 and 3, whereas the right panel shows an orthogonal 
projection from the negative y axis at the same times. 
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Polytropic scale 


K 


100 


Polytropic index 


r 


2 


Central rest-mass density 


Pc 


10-3 


Coordinate axis ratio 


rp/re 


0.7 


ADM mass 


M 


1.4906 


Rest mass 


Mo 


1.5936 


Equatorial proper radius 


Re 


12.322 


Equatorial inverse compactness 


Re/M 


7.7321 


Angular momentum 


] 


1.3192 


Normalized angular momentum 


//M2 


0.5938 


Kinetic over binding energy 


T/\W\ 


7.4792 ■ IQ-^ 


(see caption) 


Cl/ClK 


0.7536 



Table I: Parameters and integral quantities of the uniformly 
rotating polytrope used as a code test for the cubed sphere 
thirteen patches system. The quantities K, T, pc, and rp/re 
are parameters. The quantity Q is the angular velocity, while 
is the associated Keplerian velocity. Therefore, the mass- 
shedding sequence is located at fie/Cli;^ = 1 (cf. also js^ t. 



F. Uniformly rotating polytrope on the cubed sphere 
thirteen patches system 

An important equilibrium solution for general rela- 
tivistic astrophysics are rotating polytropes, since they 
provide an approximate model for relativistic stars. 
They are also particular well-suited to serve as code 
tests, since they involve a non-trivial spacetime geom- 
etry even in adapted coordinates js^, and their stability 
properties are well investigated in the case of stiff (poly- 
tropic index F = 2) and uniformly rotating solutions. 

We constructed a uniformly rotating polytropic solu- 
tion with the rns code ll60ll , assuming a stratification 
P = Kp^ with K = 100 and T = 2. The central density 
is fixed to pc = 10^^, and the ratio of polar to equatorial 
radii (in terms of the particular gauge choice in rns ll59ll ) 
is set to Tp/Ve = 0.7, which produces a rapidly rotating 
model (see Table IJl. The solution, which is represented 
on a two-dimensional grid in rns, is then mapped to ev- 
ery patch, and we apply the appropriate tensor coordi- 
nate transformations on all quantities (the four-metric, 
its derivatives, and the primitive hydrodynamical vari- 
ables). 

To specify a multi-patch system, we need to choose 
resolutions and set the free parameters tq, r\, ri (see Sec- 
tion HVEj which correspond to the location of the cube 
boundary, the spherical boundary between the seven 
patches and the six patches system, and the spherical 
outer boundary. We use tq = 1, = A and r2 = 14 
for these locations in the test case presented here, but 
we have also experimented with different values and 
obtained very similar results. The number of cells per 
patch are set to n^, where n = 20,40,80. During evo- 
lution, we do not enforce the polytropic constraint P = 
Kp^, but rather use the gamma law P = (T — 1)m. 



Quantity 20 x 20 x 20 40 x 40 x 40 80 x 80 x 80 

Central density 0.1% 0.04% 0.013% 

Rest mass 0.45% 0.16% 0.054% 

Angular momentum 1.07% 0.38% 0.13% 



Table II: Rotating neutron star on the cubed sphere system of 
thirteen patches. This table shows average errors acquired per 
dynamical time in different resolutions. 

Fig. [T7| shows the initial data mapped to the thirteen 
patches system. Note that the density outside the star 
is not exactly zero, but set to the atmosphere value (see 
Section llIB2|l since our techniques are unable to handle 
vacuum-matter interfaces. Fig. [18] shows an evolution 
of the polytrope for about 10 dynamical times^. Note 
that we use the conserved variables to show convergence: 
the three-velocity error is dominated by low-density, but 
fast material leaving the stellar surface, which does not 
carry a high amount of momentum, however. This re- 
flects the approximation we make with using an arti- 
ficial atmosphere in the first place. Therefore, conver- 
gence tests in the u' are dominated by noise. 

The evolution of the star exhibits artifact oscillations 
and a linear drift in the central density, but converges 
to the stationary exact solution. The average errors ac- 
quired by the star per dynamical time are listed in Ta- 
ble lUl 

G. Equilibrium accretion torus on a Schwarzschild 
background 

This test models a thick accretion disk in the test-fluid 
limit around a black hole. The equilibrium structure of 
these solutions has been discussed in several seminal 
papers by Fishbone and Moncrief I61|]and Abramowicz, 
Jaroszyhski, Sikora and Kozlowski ||62|, l63ll63l . For the 
purposes of this test, we use a polytropic disk with con- 
stant specific angular momentum on a Schwarzschild 
background. 

The details of how to construct these models can be 
found in the aforementioned publications, we will only 
give a quick overview here: Starting with a spacetime 
with a line element of the form 

ds^ = gttdt^ + 2gt^dtd(p + grrdr^ + geedO^ + g^dcp^, 

(19) 

we seek stationary solutions for a rotating fluid given 
by the energy-momentum tensor eqn. [T] For the four 
velocity of the fluid, we assume the form 



^ We make use of the common choice to = Re \/Re/M as a measure 
of dynamical time, where Rc is the proper equatorial circumferential 
radius and M the ADM mass of the star. 
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uf = {u\0,0,uff (20) 

and define the angular velocity G = u't' /u* and specific 
angular momentum Z = —u^/ut. From the normaliza- 
tion wf'u,, = — 1 we obtain 

{ut)-^ = g''-ig'n+gfn^. (21) 

From this quantity, and using / = const and the poly- 
tropic relation P = Kp^ , we can calculate the specific 
internal energy e = u/ p and the density Ii63il : 

e = (^i^-i^/r (22) 

Here, (Mt)o is a free parameter of the solution. The 
particular parameter values we use to construct the 
torus are K = IQ-^, F = 4/3, / = 4.5, and (Mt)o = -0.98. 
The atmospheric density is set to lO^^'^ (compared to a 
maximal density of 2 ■ 10^^). In addition, we apply a ra- 
dial coordinate transformation of the form r = exp(f), 
i.e., a logarithmic grid, to resolve the region close to the 
black hole better. The inner and outer boundaries are 
treated by extrapolation onto ghost zones as mentioned 
in SectionHnSl 

We use a cubed sphere six patches system, with 15 x 
15 X 40, 30 X 30 X 80 and 60 x 60 x 160 cells per patch 
(remember that our convention in eqn. [15] implies that 
the third local coordinate corresponds to the location 
of each sphere in the global coordinate space, whereas 
the first two coordinates roughly correspond to the an- 
gular directions on the sphere). The free parameters 
specifying the boundary location are set to rg = 6 and 
ri = 50. This places the inner radius outside of the 
horizon, which is necessary since there is a coordinate 
singularity at the horizon in these coordinates. This is 
possible since the background is not evolved and since 
there is no matter near the horizon. 

The evolution of the torus for about ten rotational 
times (in terms of the radius of maximal density) is 
displayed in Fig. |2D1 The lowest resolution case, n = 
15 X 15 X 40 per patch, is clearly under-resolved. The 
error ratio between the higher resolutions is about 0.44, 
which is only slightly better than first order convergence 
(Table |III] collects the errors for different resolutions). 
The numerical techniques we are using are "at best" 
second-order accurate, but drop to first order near max- 
ima and at the boundary to the atmosphere. In cases 
where internal turbulence is present, as is expected in 
realistic accretion disks Issl, |65|] , it will likely dominate 
the solution errors. 



Quantity 15 x 15 x 40 30 x 30 x 80 60 x 60 x 160 

Central density w 4.6% 0.055% 0.021% 
Rest mass 0.1% 0.32% 0.14% 
Angular momentum « 0.1% 0.32% 0.14% 



Table III: Equilibrium accretion torus on a Schwarzschild back- 
groimd, using the cubed sphere six patches system. This ta- 
ble shows average errors acquired per rotational time in dif- 
ferent resolutions. Note that the lowest resolution case (n = 
15 X 15 X 40 is imder-resolved. 



VI. SUMMARY 

The purpose of this paper was to demonstrate that 
multi-patch techniques can be employed to use quasi- 
spherical grids for applications in general relativistic as- 
trophysics. In addition to a standard set of techniques 
for modeling relativistic flows, we use overlapping grid 
patches with boundary ghost zones, which are syn- 
chronized to data obtained from an interpolation opera- 
tion, and a transformation associated with the transition 
map. 

The resulting scheme is able to cleanly transport 
shock fronts across patch interfaces and evolve accretion 
disks with all the advantages of a grid based on spheri- 
cal polar coordinates, but without sharing its disadvan- 
tages like artificial outer boundaries due to coordinate 
singularities, or small time steps imposed by the conver- 
gence of the latitudinal great circles towards the poles. 
In addition, spherical grids admit to choose radial and 
angular resolutions independently, which is mirrored in 
our approach, and makes it possible to use radial grids 
with sufficient radial and angular resolution near the 
black hole horizon. 

The main focus of this approach will be general rela- 
tivistic single star and star + accretion disk models, but 
there is no fundamental reason why they could not be 
employed in binary models as well. In particular, the 
inspiral phase could potentially be treated with a re- 
duced solution error compared to Cartesian mesh re- 
finement setups (see, for example, fb^). Another pos- 
sible application is to model a source with Cartesian 
mesh refinement in the domain center, but use a cubed 
sphere seven or thirteen patches system to provide outer 
boundary conditions and propagate gravitational ra- 
diation. Clearly, many of these scenarios require to 
solve Einstein's equations coupled to relativistic hydro- 
dynamics or even magnetohydrodynamics, which will 
be a topic of this series in the future. 
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Figure 13: Sod test on the six patches system. These series of 
plots demonstrate the convergence of the primitive variables 
to the solution produced by the riemann code. Each of the six 
patches has the same resolution defined by the parameter n, 
which determines the number of cells by n x n x m. 
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Figure 14: Sod test on the six patches system. The plots show a 
close-up view of the transmission of the shock front across one 
of the interfaces, using a pseudo-Schlieren plot: The function 
displayed is the norm of the density gradient, in logarithmic 
scale. The interface is between the patch boundaries indicated 
by white lines. 
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Figure 15: Sod test on the seven patches system. The left panel shows the evolution of the density for the case where each of the 
patches has resolution 80 x 80 x 80. The white lines indicate the boundaries of the seven patches which are used to cover the 
computational domain. The yellow surface is the cut of the plane z = with the grid boundaries, and the density function is 
shown at 2 = 0. The right panel shows convergence in the primitive variables. 
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Figure 16: Sod test on the thirteen patches system. The left panel shows the evolution of the density for the case where each of 
the patches has resolution 80 x 80 x 80. The white lines indicate the boundaries of the thirteen patches which are used to cover 
the computational domain. The yellow surface is the cut of the plane z = with the grid boundaries, and the density fimction is 
shown at z = 0. The right panel shows convergence in the primitive variables. 



Figure 17: Rotating neutron star on the cubed sphere system 
of thirteen patches. The initial data is produced with the rns 
code IS] then mapped to the each patch and transformed to 
the local coordinate system. These plots show the logarithm of 
the density function in the equatorial plane 2 = (left) and the 
plane x = (right). For visual clarity, the density is cut below 
10^®; the actual initial density in the atmosphere is 10^^". 
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Figure 18: Evolution of a rotating neutron star on the cubed sphere system of thirteen patches with different resolutions. The 
star is evolved up to a coordinate time of 350, which translates into about ten dynamical times (to = Re \/Re/M ^ 35.428). The 
top left panel shows the central density, the top right, middle left, and middle right panels show convergence of the conserved 
variables, and the last two plots show the errors in rest mass (bottom left panel) and angular momentum (bottom right panel). 
For a definition of these quantities see (S^ . 
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Figure 19: Equilibrium accretion torus on a Schwarzschild background, using the cubed sphere six patches system. The plots 
show, for the evolution with 60 x 60 x 160 cells per patch, the decadic logarithm of the density for coordinate times t = 0, 
1000 and 2000 (the rotational time of the torus at the locus of highest density is about 200, i.e., we are evolving for about 10 
rotational times). The left panel shows the cross-section surface x = in terms of the global, pseudo-Cartesian coordinates, and 
the right panel provides a view of the equatorial plane 2 = 0. The system is stationary, so all deviations from the initial data is 
an artifact due to the finite resolution, finite precision of floating point numbers, or the use of an artificial atmosphere (though 
small numerical errors might initiate physical instabilities in the disk). A small amount of material gets unbound and is either 
accreted into the black hole or forms a corona. In disks where a real physical effect like turbulence causes redistribution of angular 
momentum (3^, similar domains are present, though with much higher densities (note that the color map covers a range of 10 
orders of magnitude in density). 



27 



0.023 
0.022 
0.021 
0.02 
0.019 
0.018 



n = 1 5 X 1 5 X 40 
n = 30 X 30 X 80 
n = 60 X 60 X 1 60 
reference (see caption) 



100 
10 

S 1 
Q 

0.1 

0.01 



n = 1 5 X 1 5 X 40 
n = 30 X 30 X 80 
n = 60 X 60 X 1 60 




1e-04 

0.03 r 

0.02 - 

0.01 - 


-0.01 - 

-0.02 - 

-0.03 - 

-0.04 - 




500 



1000 

t 



1500 



0.01 

2000 



500 



1000 

t 



1500 



2000 



n = 1 5 X 1 5 X 40 
n = 30 X 30 X 80 
n = 60 X 60 X 1 60 
reference (see caption) 




500 



1000 
t 



1500 



2000 



0.02 - 

0.01 - 

I 

-0.01 - 

-0.02 - 

-0.03 - 

-0.04 - 




n = 1 5 X 1 5 X 40 
n = 30 X 30 X 80 
n = 60 X 60 X 1 60 
reference (see caption) 




500 



1000 
t 



1500 



2000 



Figure 20: Equilibrium accretion torus on a Schwarzschild background, using the cubed sphere six patches system. The disk is 
evolved up to a coordinate time of 2000, which translates into about ten rotational times at the location of maximal density. The 
top left panel shows the maximal density, the top right, middle left, and middle right panels show convergence of the conserved 
variables, and the last two plots show the errors in rest mass (bottom left panel) and angular momentum (bottom right panel). 
For a definition of these quantities see [59]. (The "reference" solution mentioned in the text is the value corresponding to the 
initial data on grid of 120 x 120 x 320 cells per patch.) Note that the case n = 15 x 15 x 40 is clearly under-resolved. 



